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We calculate the angular two-point correlation function of ultra-high energy cosmic rays (UHECR) ob- 
served by AGASA and Yakutsk experiments. In both data sets, there is a strong signal at highest energies, 
which is concentrated in the first bin of the size of the angular resolution of the experiment. For the uniform 
distribution of sources, the probability of a chance clustering is 4 x 10"". Correlations are absent or not 
significant at larger angles. This favors the models with compact sources of UHECR. 
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1. The measurements of the flux of UHECR at ener- 
gies of order 10^° eV provide compeUing evidence 
,of the absence of the Greisen-Zatsepin-Kuzmin (GZK) 
cutoff H]. The resolution of this puzzle seems to be 
impossible without invoking new physics or extreme as- 
trophysics. All models suggested so far can be classified 
in three groups, according to the way the GZK cutoff 
is avoided: i) "nearby source" ii) weak interaction with 
CMB iii) bump in the injection spectrum. 

The possibility i) assumes that a substantial fraction 
of the observed UHECR comes from a relatively nearby 
'source(s) and thus is not subject to the GZK cutoff. 
,This idea may be realized in different ways, examples 
being the models of decaying superheavy dark matter 
or models in which UHECR emitted by nearby source(s) 
propagate diffusively in the galactic [Q or extragalactic 
1^ magnetic fields. Although models of this type gener- 
ically predict large-scale anisotropy j|, D, they might 
still work. 

In the option ii) the GZK cutoff is eliminated (or 
shifted to higher energies) by assuming weak or non- 
standard interaction of primary particles with the cos- 
mic microwave background. This may happen, for in- 
stance, if primary particles are neutrinos , hypotheti- 
cal light SUSY hadrons ||] , or if the Lorentz invariance 
is violated at high energies §. The possibility iii) can 
be realized in models which involve topological defects 
||lo| or in some models where primary particles are neu- 
trinos 

Existing data hint also at another important fea- 
ture of UHECR, namely, the clustering at small an- 
gles . The AGASA collaboration has reported three 
doublets and one triplet out of 47 events with energies 
E > 4 X 10^^ eV, with chance probability of less than 



1% in the case of the isotropic distribution The 
world data set has also been analyzed; 6 doublets and 2 
triplets out of 92 events with energies E > A x 10^^ eV 
were found [|^, with the chance probability less than 
1%. 

If not a statistical fluctuation, what does the clus- 
tering imply for models of UHECR? There are two pos- 
sible situations: either clustering is due to the existence 
of point-like sources, or it is a result of variations in the 
flux of UHECR over the sky (the regions of higher flux 
are more likely to produce clusters of events |l^). In 
the first case, the models which involve the diffuse prop- 
agation of UHECR are excluded. This case also implies 
that there is no defocusing of UHECR in the magnetic 
fields during their propagation. Such defocusing occurs 
even in a regular (e.g., galactic) magnetic field, since dif- 
ferent events in a cluster have different energies. Thus, 
one can put bounds on the charge of the primary parti- 
cles. For the extra-galactic rays, knowing that primary 
particles are charged would imply direct bounds on the 
extra-galactic magnetic fields. 

In the second case, the regions of higher flux may 
reflect higher density of sources as in the models of su- 
perheavy dark matter where they would correspond to 
dark matter clumps in the halo. Alternatively, they may 
be due to the effects of propagation such as defocusing 
in magnetic flelds or magnetic lensing |p^ . 

In order to determine which of these two cases fits 
the present experimental data better, it is not enough to 
know the probability to have a certain number of clus- 
ters. In this respect previous analyses ^ are 
not sufficient. One has to find the angular correlation 
function. This is the approach we accept in this paper. 
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2. The two-point correlation function for a given set 
of events is defined as follows. For each event, we di- 
vide the sphere into concentric rings (bins) with fixed 
angular size (say, the angular resolution of the experi- 
ment). We count the number of events falling into each 
bin, sum over all events and divide by 2 to avoid dou- 
ble counting, thus obtaining the numbers A^^. We repeat 
the same procedure for a large number (typically 10^) of 
randomly generated sets and calculate the mean Monte- 
Carlo value Nf^'^ and the variance cr^^ for each bin in a 
standard way. The correlation function can be defined 
as fi = Ni/Nf^'^ — 1. A deviation of fi from zero in- 
dicates the presence of the correlations on the angular 
scale corresponding to i-th bin. 

The correlation function fi fluctuates. In order to 
see whether its deviation from zero is statistically signif- 
icant, we define the ratio ri = {Ni - Nf^'^)/af^^, which 
shows the excess in the correlation function as compared 
to the random distribution in the units of the variance. 
With enough statistics, this quantity becomes a good 
measure of the probability of the corresponding fiuctu- 
ation. 

The Monte-Carlo events are generated in the hori- 
zon reference frame with the geometrical acceptance 
dn (X cos 9 z sin 6 zddz, where 9z is the zenith angle. Co- 
ordinates of the events are then transformed into the 
equatorial frame assuming random arrival time. We 
restrict our analysis to the events with zenith angles 
9z < 45° for which the experimental resolution of ar- 



rival directions is the best 1 14 



If clusters at highest energies are not a statistical 
fiuctuation, one should expect that the spectrum con- 
sists of two components, the clustered component tak- 
ing over the uniform one at a certain energy. The cut 
at an energy at which the clustered component starts to 
dominate should give the most significant signal. Mo- 
tivated by these arguments, we calculated the proba- 
bility of chance clustering as a function of energy cut. 
We present here the results for the AGASA ||l^ and 
Yakutsk data sets (other experiments are discussed 
in Section 3). For these simulations we took the bin size 
equal to 2.5° and 4° for AGASA and Yakutsk, respec- 
tively, which is the quoted (see e.g. [^, |lj|) angular 
resolution of each experiment multiplied by -\/2. The 
results are summarized in Fig. ^ which shows the prob- 
ability to reproduce or exceed the observed count in the 
first bin, as a function of the energy cut. AGASA curve 
starts at = 4 X 10^^ eV because the data at smaller 
energies are not yet available. Yakutsk has much lower 
statistics. Both curves rapidly rise to 1 in a similar way 
when the statistics becomes poor. They suggest that 
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FIG. 1. Probability to match or exceed the observed 
count in the first bin as a function of energy for the 
random distribution of arrival directions. 



the optimum energy cut is higher than can be imposed 
at present statistics. 

The difference between our results and those of 
Ref. Q (cf. Fig. of this paper and Fig. 12 of Ref. |l|) 
is due to two reasons. First, 10 more events with 
E > 4 X 10^^ have been observed |18 which bring a 
new doublet. Second, and more important, we calcu- 
late a different probability. The difference arises when 
there is a triplet or higher multiplets in the data. In 



our approach a triplet is equivalent to 3 or 2 doublets, 
depending on the relative position of the events (com- 
pact or aligned), while higher multiplicity clusters effec- 
tively have larger "weight". In Ref. the probabil- 
ities of doublets and triplets are calculated separately; 
the probability of doublets is defined in such a way that 
a triplet contributes as 3/2 of a doublet. The drawback 
of this method is that the probabilities of doublets and 
triplets are not independent, and it is not clear how to 
combine them. Triplets and higher multiplicity clusters 
are better accounted for in our method, and the proba- 
bility of chance clustering which we get is lower than in 
Ref. |l|. 

Correlation functions calculated with the energy 
cuts corresponding to the lowest chance probability is 
shown in Fig. |^. Both AGASA and Yakutsk correla- 
tion functions have substantial excess in the first bin. 
The peak in AGASA curve corresponds to 6 doublets 
(of which 3 actually form a triplet) out of 39 events. 
The peak in Yakutsk curve corresponds to 8 doublets 
(of which 3 also come from a triplet) out of 26 events. 

Since the number of events in the first bin is not 
large, its distribution is not well approximated by the 
Gaussian one, and the deviation in the units of the 
variance is not a good measure for the probability of 
fluctuations. We calculated the probability directly by 
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FIG. 2. Angular correlation functions of UHECR with 
binning angles and cuts in energy quoted in the text. 



counting, in the Monte-Carlo simulation, the number of 
occurrences with the same or larger number of events in 
the first bin. Probabilities in the minima are small (see 
Fig. |l|), so we have recalculated them with 10^ Monte- 
Carlo sets. 

As the lowest probabilities were obtained by scan- 
ning over the energy, one may argue that they have to 
be multiplied by the number of steps in the scan. This 
is not, however, correct because the results at different 
energy cuts are not independent: higher energy set is 
a subset of the lower energy one. As can be seen from 
Fig. 1^, the chance probability for AGASA is lower than 
10~3 in the whole energy range (4 — 5) x 10*^ eV, re- 
gardless of the number of steps in the scan. There still 
may be a correction factor. To estimate it we made the 
following numerical experiment. For 10^ randomly gen- 
erated sets of events we have performed exactly the same 
procedure as for the real data, i.e. scanned over energies 
and obtained 10^ different minimum probabilities. We 
found that the probability less than 10~^ occurred 27 
times, while the probability less than 10"^ occurred 3 
times. Thus we conclude that the correction factor is 
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of order 3. This factor is included in the final results 
which are presented in Table. 1. 

We now turn to the determination of the angular 
size of the sources. To this end we calculate the depen- 
dence of the probability to have the observed (or larger) 
number of events in the first bin on the bin size. This 
dependence is plotted in Fig. ^. Jumps in the curves 
occur when a new doublet enters the first bin. Despite 
fluctuations, one can see that the minimum probabil- 
ity corresponds roughly to 2.5° and 4° for AGASA and 
Yakutsk, respectively. These numbers coincide with the 
angular resolutions of the experiments, as is expected for 
sources with the angular size smaller than the experi- 
mental resolution. Remarkably, there are no doublets 
in the AGASA set with separations between 2.5° and 
5° , while for the the extended source of the uniform lu- 
minosity one would expect 4 times more events within 
5° as there are within 2.5°. Thus, we conclude that the 
data favor compact sources with angular size less than 
2.5°. 

If primary particles are charged, actual positions of 
sources differ from the measured arrival directions be- 
cause of the deflection in the Galactic magnetic field 
(GMF). If the clustering is attributed to real sources, it 
should not disappear but improve when the correction 
for GMF is taken into account. We have simulated the 
effect of such correction making use of the GMF mod- 
els summarized in Ref. [Q. For the charge Z = 1 and 
BSS_A model the peak in the first bin does not change 
significantly; one cannot discriminate between this case 
and the case of neutral particles. The peak becomes 
small at Z = 2, and disappears at larger Z for all GMF 
models of Ref. [01. 



3. The other two UHECR experiments, Haverah Park 
(HP) and Volcano Ranch (VR), do not see significant 
clustering [|||. With the energy cut £: > 2.4 x lO^^ eV 
and the bin size 4°, the HP data contain 2 doublets at 
1.8 expected, while VR data contain 1 doublet at 0.1 ex- 
pected with isotropic distribution. Let us estimate the 
combined probability of clustering in all experiments as- 
suming independent Poisson distributions. The number 
of observed doublets in AGASA and Yakutsk data are 
6 and 8, respectively, while 0.87 and 2.2 are expected 
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FIG. 3. Probability to have observed count in the first 
bin as a function of the bin size. Cuts in energy corre- 
spond to minima of Fig. 1. 

Table 2 





Ntot 


observed 


expected 


probability 


AGASA 


39 


6 


5.4 + 0.6 




Yakutsk 


26 


8 


2.9 + 1.6 


0.09 


HP 


32 


2 


4.0+1.8 


0.07 


VR 


10 


1 


0.7 + 0.1 


0.55 



(these "effective" expected numbers of doublets are cal- 
culated from the condition that probabilities of Table 1 
are reproduced, i.e. "penalty" for the energy scan is 
included). Thus, 17 doublets are observed at 4.97 ex- 
pected, which corresponds to the Poisson probability 
2 X 10^'''. If HP data are excluded, the probability be- 
comes 1 X 10~^, while with both HP and VR data ex- 
cluded the probability is 4 x 10~^. 

It is extremely unlikely that the clustering observed 
by AGASA and Yakutsk experiments is a result of a 
random fluctuation in an isotropic distribution. Rather, 
the working hypothesis should be the existence of some 
number of compact sources which produce the observed 
multiplets. Is this hypothesis consistent with HP and 
VR data? For a given experiment, the expected num- 
ber of clusters is determined by the total number of 
events poj (see also ref.|^); at small clustering it scales 
like A^toT |§. Taking AGASA data as a reference (6 
doublets observed, 5.4 expected from sources and 0.6 
expected from chance clustering) allows to estimate the 
expected number of doublets in other experiments by 
adding the doublets expected from sources and the dou- 
blets expected from the uniform background (calculated 
in the Monte-Carlo simulation). The results are sum- 
marized in Table 2, together with corresponding Poisson 
probabilities. 



FIG. 4. Observed clusters in Galactic coordinates. 

All experiments are roughly consistent with the as- 
sumption that the number of sources is such that they 
produce 5.4 doublets out of 39 events in average. Note 
that if HP data are discarded the agreement be- 
tween other experiments can be made better. 

According to our simulations, the mean numbers 
of chance doublets are 0.6 and 1.6 for AGASA and 
Yakutsk, respectively. Therefore, most of the clusters 
in AGASA and Yakutsk data are likely to be due to 
real sources. In Fig. |4| we plot these clusters in the 
Galactic coordinates (small and large circles correspond 
to AGASA and Yakutsk events, respectively). Positions 
of triplets arc indicated by arrows. The set of AGASA 
events with £' > 4 x 10^^ eV and Yakutsk events with 
E > 2.4 X 10^^ eV is a suitable choice for the search of 
correlations with astrophysical objects. 

To summarize, the clustering of UHECR is statisti- 
cally significant and favors compact sources. This places 
further constraints on models which can resolve the puz- 
zle of the GZK cutoff. Those models which involve large 
extragalactic magnetic fields, Ref. as well as models 
with heavy nuclei as primaries, e.g. are disfavored 
because they assume total isotropisation of original ar- 
rival directions of UHECR. If violation of the Lorentz 
invariance is the solution of the GZK puzzle, and pri- 
maries are protons, our results place extremely strong 
limit on the extragalactic magnetic fields. Regarding 
the models of decaying superheavy dark matter, it is 
important to calculate |2^ the angular correlation func- 
tion predicted by these models and compare it to Fig.|| 
in order to see if the clumping on subgalactic scales can 
be responsible for the clustering of UHECR. 
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